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A PERTURBATIVE ANALYSIS OF QUASI-RADIAL DENSITY WAVES IN 

GALACTIC DISKS. 
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RESUMEN 
ABSTRACT 

The theoretical understanding of density waves in disk galaxies starts from 
the classical WKB perturbative analysis of tight-winding perturbations, the key 
assumption being that the potential due to the density wave is approximately ra- 
dial. The above has served as a valuable guide in aiding the understanding of both 
simulated and observed galaxies, in spite of a number of caveats being present. The 
observed spiral or bar patterns in real galaxies are frequently only marginally consis- 
tent with the tight-winding assumption, often in fact, outright inconsistent. Here 
we derive a complementary formulation to the problem, by treating quasi-radial 
density waves under simplified assumptions in the linear regime. We assume that 
the potential due to the density wave is approximately tangential, and derive the 
corresponding dispersion relation of the problem. We obtain an instability criterion 
for the onset of quasi-radial density waves, which allows a clear understanding of the 
increased stability of the higher order modes, which appear at progressively larger 
radii, as often seen in real galaxies. The theory naturally yields a range of pattern 
speeds for these arms which appears constrained by the condition f2p < JIq ± n/m. 
For the central regions of galaxies where solid body rotation curves might apply, 
we find weak bars in the oscillatory regime with various pattern speeds, including 
counter rotating ones, and a prediction for Q.p to increase towards the centre, as 
seen in the rapidly rotating bars within bars of some numerical simulations. We 
complement this study with detailed numerical simulations of galactic disks and 
careful Fourier analysis of the emergent perturbations, which support the theory 
presented. 

Key Words: INSTABILITIES STELLAR DYNAMICS GALAXIES: 
KINEMATICS AND DYNAMICS GALAXIES: STRUC- 
TURE GALAXIES: SPIRAL 



1. INTRODUCTION 

The development of short-wavelength, tight- 
winding disk wave dynamics, notably by Lin & Shu 
(1964), Shu (1970) and Toomre (1969), still forms 
the basic analytical substrate used to understand 
galactic structure observations and numerical sim- 
ulations of galactic dynamics. This WKB perturba- 
tive analysis exploits the fact that for tight-winding 
spirals, the long range character of the gravitational 
potential due to the density wave one is modeling 
disappears, the response becomes local, and an ana- 
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lytical formulation to the problem becomes feasible. 
The dispersion relation that results has been used 
successfully over the decades to gain some under- 
standing of the expected physical scalings between 
pattern speed, the orbital frequency, the epicycle fre- 
quency and the mass density of galactic disks. 

However, several weak points in the application 
of the classical WKB density wave dynamics to real 
galaxies have been recognised since they where orig- 
inally proposed. First, one would like spiral arms 
to be very tightly wound, making the approximation 
on the radial dominance of the perturbed potential 
clearly valid. For typical early type spirals, spiral 
arms are clearly tightly wound, but often one finds 
late type spirals with almost radial 'spiral' arms. 
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The spiral patterns of many galaxies appear to be 
on the limit of applicability for this approximation, 
and many, galactic bars included, are outright out of 
it. 

One of the strong points of the density wave the- 
ory was the substantial easing of the winding prob- 
lem, which can actually be used to rule out any mate- 
rial arms interpretation for spiral structure in galax- 
ies, provided one seeks a long lived structure. 

Here we explore a complementary, and in a cer- 
tain sense, perpendicular approach. We shall assume 
that the density wave being treated is quasi-radial, 
and that it represents only a small disturbance on 
the overall potential. Under this assumptions, the 
potential of the disturbance will only be felt near 
the density enhancement, and will be tangential, di- 
rected towards the density enhancement itself, par- 
ticularly if one assumes no strong radial gradients in 
the amplitude of the density wave itself. Treating 
only the circular motions of a coUisionless compo- 
nent, and ignoring the vertical gradients of all quan- 
tities involved, lands us in the thin disk approxima- 
tion for stellar disks, a scenario often used in galactic 
dynamical works and in studies of galactic density 
waves in analytical developments, orbital structure 
analysis or N-body simulations, being Hernandez & 
Cervantcs-Sodi (2006), JalaU & Hunter (2005), Ma- 
ciejewski & Athanassoula (2008) and O'Neil & Du- 
binski (2003) examples of the above. 

This allows an approximate analytical treatment 
of the problem which we present here in section (2). 
A dispersion relation for the problem appears, the 
consequences of which we explore under two ide- 
alised galactic regimes, flat rotation curve, and solid 
body rotation, with the aim of sampling the resulting 
physics under the approximate dynamical conditions 
where one finds spiral arms and bars. As in the case 
of the Lin-Shu formalism, we begin by treating only a 
stellar disk, thinking of modeling the old stellar pop- 
ulation of a present day galaxy. The development is 
idealised and remains within the linear regime, for 
quasi-radial density waves, for a disk made up of 
close to circular orbits. While lacking the detail of 
more advanced studies (e.g. Evans & Read 1998, 
Polyachenko 2004, Jalali & Hunter 2005), it has the 
advantage of identifying clearly a physical instability 
criterion which captures some of the details missed 
by the Lin-Shu formalism, highlighting the role of 
stellar diffusion. As an independent confirmation of 
the ideas presented, we also run an extensive grid 
of numerical models, and analyse the details of the 
small amplitude density perturbations which appear. 
It is encouraging that the results of these direct nu- 



merical simulations validate the general results of the 
ideas developed in this study. 

In section (3), for the case of the fiat rotation 
curve regime, we derive a dispersion relation for 
the problem, and a corresponding stability criterion 
for the onset of quasi-radial density perturbations, 
which is seen to be clearly satisfied for any reason- 
able galactic disk, for low values of m and Q, the 
symmetry number for the pattern and Toomre's pa- 
rameter. The case of a solid body rotation regime 
is also treated, where we show that for weak bars 
in the oscillatory regime of the dispersion relation, 
relatively constant values for fJp appear quite natu- 
rally, at values which represent 'slow', 'fast' or even 
counter-rotating bars, for natural values of m and 
Q. The extent of these bars being naturally limited 
by the corotation condition. A detailed numerical 
model is presented in section (4), including a careful 
Fourier analysis of the density perturbations which 
appear, yielding results in consistency with the the- 
ory developed. Section (5) presents our conclusions. 

2. PHYSICAL SETUP 

In this section we develop the simplified physical 
model for quasi-radial density waves. As opposed 
to the Lin-Shu formalism, designed to model tight- 
winding waves, we assume that the derivatives of the 
potential due to the perturbation are dominated by 
the angular component, rather than the radial one. 
We shall start from the angular moments equation 
written in cylindrical coordinates {R,0,z), in a ref- 
erence frame which rotates with a constant angular 
frequency flp parallel to the z axis, but without as- 
suming azimuthal symmetry. In spite of the multi- 
component nature of a galactic system, we shall be- 
gin by ignoring the gaseous component, in an at- 
tempt to model the imderlying density perturba- 
tion represented by loosely-wound arms in the old 
stellar population. We shall also assume no ver- 
tical or radial velocities, in modelling weak radial 
arms which only result in weak angular distortions 
in both the density and velocity fields. Taking also 
an isotropic velocity ellipsoid with no vertex devia- 
tion, i.e. aij = Sija, results in the following equation 
(e.g. Vorobyov & Theis 2006): 

dt R 09 09 09 ^ ' 

In the above equation <^^f = <I> — x r|^/2 and 
S is the mass surface density of the disk. Taking 
the quasi-radial axms we are interested in modelling 
as a perturbation on the axisymmetric centrifugal 



QUASI-RADIAL DENSITY WAVES IN GALAXIES 



3 



equilibrium solution gives rise to the following set of 
conditions at every radius: 

S = So + eSi(^) 

ve = vo{R) + €Vi{0) 
a = ao + ecTi (6) 

In the above $o is the axially symmetric gravita- 
tional potential, which cancels the centrifugal force 
fixing v(r) = (0, Rifto-Qp), 0), with QoiR) the cen- 
trifugal equilibrium orbital frequency $1 will hence 
be the potential perturbation due to the quasi-radial 
arms, resulting in angular perturbations Si(^), vi{6) 
and (Ji[0.) The form taken for the velocity implies 
no further inertial terms are necessary. The unper- 
turbed state trivially satisfies eq.(l) for the standard 
axisymmetric centrifugal equilibrium state. 

To first order in the perturbation, eq.(l) now be- 
comes: 



REo^ + R^An^ 

+{R^A^n + 4)^ 



(2) 



In the above equation we have introduced Af2 = 
(ilo — fip), and have omitted making explicit the ra- 
dial dependence of all the above variables, which re- 
mains implicit. 

We now introduce an isothermal equation of state 
(e.g. Binney & Trcmaine 1987) or equivalently the 
conservation of phase-space density E/a^ = cte. to 
eliminate ai for Si through: 



(5) 



We will change the dependencies of the above 
equation on vi for dependencies on Ei through the 
use of the continuity equation, which reads: 



^E 1 d{Ri:vn) 1 d{J:ve) _ 

dt R dR R de ■ ^ ' 

Introducing the perturbation conditions with 
only tangential velocities into the above continuity 
equation yields: 

dt R de de ^ ' 

From this last equation we can solve for dvi/d6 
to construct: 



d'^vi RAn a^Ei R S^Ei 



Eo 9612 
RAn a2Ei 



Eo dedt 

R a2Ei 



(8) 



(9) 



dedt Eo dedt Eo dt^ 

The last two above are now used to replace the 
dependence on vi appearing in eq.(5) for a depen- 
dence on El, giving: 



-R 



2 a^Si 

-^O-ggT- 



2R^AVl 



aedt 



+ ( 



^0-2 



i?2A2n) ^ 



(10) 

Finally, we relate Ei and $1 through the first 
order Poisson equation of the problem, which under 
the assumption of tangential forces due to the almost 
radial perturbation dominating over the other two 
directions gives: 



with which eq.(2) becomes: 



(3) 



1 92*1 



i?2 de^ 

allows to write eq.(lO) as: 



(11) 



i?Eof^ + i?2AOf^ 



{R^A^n - 



3"oj de 



2RAnT,o-. 



(4) 



We sec the assumption of eq.(3) in the numerical 
constant of the fifth term in the above equation. Any 
other similar equation of state would only change 
this constant slightly, leaving all conclusions quali- 
tatively unchanged, and only modified to within a 
small numerical factor of order unity. The partial 
angular derivative of the above equation reads: 



-R 



.2 

47rGpoSi 



R^A^^) 



(12) 

In the above equation we have used h, the height 
scale of the disk, as h = Eo/po = Ei/pi. 

We now turn to well known galactic structure 
scalings to re- write the right hand side of eq.(12) in 
terms of more helpful dynamical variables, specifi- 
cally, we shall make use of vertical virial equilibrium 
in the disk and Toomre's stability criterion. 
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ttGSo 



= Q. 



(13) 



In the above = 4f2§ + 29.oR{d9.o / dR) is the 
epicycle frequency as functions of the radius. The 
star formation processes in the disks of spiral galax- 
ies have often been thought of as constituting a self 
regulated cycle. The gas in regions where Q < a cer- 
tain critical vahic finds itself in a regime where lo- 
cal self-gravity, appearing through (GSo) in cq.(13), 
dominates over the combined effects of rotational 
shears (or equivalently global tidal forces) and the 
stabilising thermal 'pressure' effects of the product 
(kct) . This leads to collapse and the triggering of star 
formation processes, which in turn result in signifi- 
cantly energetic events. The above include radiative 
heating, the propagation of ionisation fronts, shock 
waves and in general an efficient turbulent heating 
of the gas media, raising a locally to values result- 
ing in Q > a certain critical threshold. This restores 
the gravitational stability of the disk and shuts off 
the star formation processes. On timescales longer 
than the few x 10 million years of massive stellar life- 
times, an equilibrium is expected where star forma- 
tion proceeds at a rate equal to that of gas turbulent 
dissipation, at time averaged values of Q Qcrit- 
Examples of the above can be found in Dopita & 
Ryder (1994), Koeppen et al. (1995), Firmani et al. 
(1996) and Silk (2001). 

The preceding argument applies to the gaseous 
component, not the stellar one, which is the one we 
are interested in modeling, particularly the old stel- 
lar component. However, if stars retain essentially 
the velocity dispersion values of the gas from which 
they formed, one also expects Q ^ Qcrit for the stars. 
Alternatively, if we consider the dynamical heating 
of the stellar populations in the disk through encoun- 
ters with giant molecular clouds or the spiral arms 
themselves, slightly larger values of Q would be ex- 
pected for the stars than for the gas. We shall there- 
fore take eq.(13) as valid throughout the disk (e.g. 
O'Neill & Dubinski 2003, Maciejewski & Athanas- 
soula 2008), without specifying any particular val- 
ues of Q at this point. Imposing viral equilibrium 
in the vertical direction in the disk (e.g. Binney & 
Tremaine 1987) yields: 



h = 



(14) 



Since h = Yip/ pQ, we can replace the dependence 
on So for one on po in eqs.(13) and (14), and solve 
for {(Jo/h) in both eqs.(13) and (14) to obtain: 



(15) 



Equation (15) serves to re-write the right hand 
side of eq.(12) and obtain: 



2i?2A17^ + (|ag-i?2A2^l) 



902 



= -8i:r(e)' 



(16) 

Equation (16) is now a second order partial dif- 
ferential equation for the temporal and angular vari- 
ations of the density of the perturbation. This is 
valid at each radius, and shows dependence on both 
the gravitational dynamics of the system through k 
and ri, and on the local stability and structure of 
the disk through Q and uo. Now we impose periodic 
sohitions in 6 for the perturbation in density, with 
an oscillatory or growing time dependence: 



El oc e 



-i{m9—u)t) 



(17) 



with m an integer, which when introduced into 
eq.(16) gives: 



(18) 



Equation (18) is the dispersion relation for the 

problem. For a density wave of the type being de- 
scribed to develop, we require for lu to be imaginary, 
i.e., 



2k /5 



1/2 



mcTo 



(19) 



This last is an instability criterion which gives the 

values of ctq below which a certain m-symmetry pat- 
tern will begin to grow, for a given rotation curve and 
value of Q. We can also obtain directly the ranges of 
values J7p is expected to take, considering that this 
limit values will be given by the condition w = in 
equation (18). Since the pattern will develop only 
when the second term in the R.H.S. of eq.(18) domi- 
nates over the first, we can start by ignoring the first 
term on the R.H.S. of this equation. For example, in 
a flat rotation curve region, with a circular velocity 
of 220 km/ s and a cr = 30km/s, even with a Q = 3, 
the first term above is only 1.7 % of the second for 
m = 1, one has to go to m = 4 for this first term 
to pass 20 % of the second. We can hence estimate 
the limit values of Clp by neglecting this first term to 
give: 
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= ±2^2 



Q' 



(20) 



Remembering that AQ = Qq — Q^, we obtain: 



'2^2^ 



(21) 



Since Q is of the order of 2\/2, the above equa- 
tion is a direct analytical explanation for the ranges 
of pattern speeds seen in numerical experiments, in- 
dependent of the swing amplification hypothesis. 

Going back to equation(19), we sec that the lower 
values of m are easier to excite than the higher ones, 
as often found in numerical simulations, providing an 
explanation for observed galaxies with strong, high 
pitch angle arms, being preferentially m = 2 sys- 
tems. Indeed, Martos et al. (2004) showed that the 
TO = 4 optical spiral pattern in the Milky Way can 
form as a consistent hydrodynamical response of the 
gas to an underlying m = 2 pattern in the old stel- 
lar population, as traced by COBE-DIRBE K-band 
photometry. 

Remembering that k/Q is really just a place- 
keeper for {GpoY^"^, the instability criterion of 
eq.(19) in the original terms reads: 



R 



> 



(- 



1/2 



1 

Gpo 



1/2 



(22) 



In this terms, the instability criterion of eq.(22) 
is just a comparison between the local gravitational 
timescale of the disk, the free-fall timescale, and the 
inter-arm diffusive crossing time. That is, the insta- 
bility will not develop if diffusion is such that the 
perturbation is blurred out before it can condense 
gravitationally. The lower the number of arms, the 
easier it is to excite the pattern, as in high m pat- 
terns, it is easier for stars to diffuse from one arm to 
another, blurring and erasing the proposed pattern. 
To conclude, we identify the dimensionless number 



/Sy/^mQg / 5 V'^ ma 
\24j kR \12-kGp) R ^ ' 

as the critical indicator for an n-armed quasi ra- 
dial density wave to develop, with the instability 
threshold appearing at P < 1. 

The consequences of eq.(23) for spiral arms and 
bars will be explored in the following section, under 
two idealised regimes, flat rotation curve, and solid 
body rotation, for natural values of the parameters 
m and Q. 



1/2 



Pm.= \ — 



3. CONSEQUENCES FOR GALACTIC DISKS 

In this section we shall trace some of the con- 
sequences for the onset of the instability presented, 
in the context of thin galactic discs. For simplicity, 
in this first treatment, we shall consider only two 
distinct regimes for the rotation curve of a galaxy, 
a strictly flat rotation curve regime, of relevance to 
the development of galactic arms, and a strictly solid 
body rotation curve region, appropriate for the mod- 
eling of dynamics of bars. The above with the inten- 
tion of developing as clear a physical understanding 
of the effect being introduced as possible. A more 
detailed treatment of the problem using more real- 
istic rotation curves and detailed numerical galactic 
models, appears in section (4). 

3.1. Galactic Arms in Flat Rotation Curves 

We start this section with a consideration of the 
flat rotation curve limit of the relations derived pre- 
viously, as a suitable first approximation to the dy- 
namics of a galaxy in the region over which galac- 
tic arms are seen. In a flat rotation curve regime, 
K = 2^/^f2, and therefore, eq.(19) becomes: 



->(5 

(To V3 



1/2 



nQ 



(24) 



In real galactic systems, typically with 5 < 
V^/(To < 8, we sec that in a flat rotation curve regime 
the instability criterion will always be satisfled for 
TO = 2, for example, inaQ = lor(3 = 2 disk. For 
TO = 4 it is still easy to satisfy the instability crite- 
rion for the onset of a strong perturbation, for Q <2>, 
but higher values of m would be hard to accommo- 
date. Indeed, in the detailed stability analysis of 
simple power-law disks in idealised rotation curves 
of Evans & Read (1998), it is shown explicitly that 
lower TO modes are always the most unstable ones. 
Actually, from the radial dependence of Pmi it can be 
seen that higher to patterns will be easier to excite 
as one moves further out in a galactic disk, in con- 
sistency with the bifurcations often seen in galactic 
spiral arms, patters which tend to be characterised 
by higher order symmetries as one moves to larger 
radii. 

It is interesting to note that since Q scales lin- 
early with (T, P„j which scales with Qc, is a quadratic 
function of a. This last means that for small values 
of (T, Pm < Q, with the opposite being true for large 
values of a. By setting P^ = Q we obtain 



ma 



1/2 



(25) 
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for low values of m this condition will result in a 
value of a larger than that present in normal galac- 
tic disks, and hence we see that typically galaxies 
will lie in the region P„j < Q. This is interesting, 
as perhaps the identification of values of Q '^^^ 3 as 
necessary to stabilise galactic disks reported for nu- 
merical experiments, might correspond to cases of Q 
stable, P unstable disks, with quasi-radial density 
waves arising, and then being wound up into tight 
spirals. 

Remembering again that (k/Q)^ is a proxy for 
Gpo7 we can compare equation (18) to the classi- 
cal WKB equivalent relation, e.g. eq.(6.61) Binney 
& Tremaine (2008). We see that the term includ- 
ing (Jo in equation (18) appears in replacement of 
the standard k term. By considering a purely radial 
perturbation with a purely tangential force field and 
resulting matter flows, we have lost the contribu- 
tion of differential rotation, in what is essentially an 
analysis at constant radius, in favour of a dynamical 
pressure term through the stellar velocity dispersion, 
which docs not appear in the classical criterion. 

3.2. Central Bars in Solid Body Rotation Curves 

In moving to central galactic bars we now go to 
a solid body rotation regime. The above considera- 
tion leads to the use of to = 2 in what follows. In 
this regime, k = 2^0, with fi a constant, and the 
instability criterion of eq(19) is modified by the sub- 
stitution of a 6 for the 3 in the denominator of the 
square root term in the right. In terms of flo this 
reads: 



mo 



> 



1/2 



(26) 



This time we see that the criterion for the on- 
set of a quasi-radial perturbation is guaranteed to 
fail inwards of a certain radius, as the term on the 
left scales linearly with radius. Perturbations in the 
above regime will hence be in the weak oscillatory 
linear regime of eq.(18). The limit values for ftp for 
these weak pulsating bars can be found from eq.(18) 
setting uj — 0, which under solid body rotation gives: 



^ = 1± 

n 



3 \RnJ \Q 



-,1/2 



(27) 



We find a slow and a fast solution, as a func- 
tion of the radial coordinate, fig, ctq and Q. Fig- 
ure (1) now gives values of i}p/Q,o, as a function of 
the dimensionless radial coordinate v — RVLq/c^), for 
3 values of Q, 1, 2 and 3, inner, middle and outer 



0.2 



Fig. 1. Values of Q,p/Q. for stationary weak bars, as 
a function of the dimensionless radial coordinate v = 
R^l/ao, for Q=l, inner curves, Q=2, middle curves and 
Q=3, outer curves. The vertex of the curves defines the 
extent of the weak bars in question, always at the coro- 
tation radius. 

curves, respectively. Notice that for large values of 
V, relatively constant values of Qp appear, for val- 
ues of Q > 2. In going towards the centre however, 
the divergence towards R = becomes evident, and 
large variations in Vlp appear. This last feature might 
explain the appearance of bars within bars, which 
exhibit increasingly rapid pattern rotation rates in 
going towards the central regions (e.g. Heller et al. 
2007, Maciejewski & Athanassoula 2008 and Shen 
& Debattista 2008). We see that the lower branch 
of the curves, corresponding to the minus sign in 
eq.(28), can have negative values over much of its 
extent. This corresponds to counter-rotating bars, 
where the stars are not counter rotating themselves, 
but the density enhancement is. 

Also, we see that replacing the inequality for an 
equality in eq.(27) we can solve for Rm, the max- 
imum extent of a weak bar of the type being de- 
scribed. This last is given by the vertex of the 
curves in figure(l), and as can be seen from eq.(27), 
will scale linearly with Q. From eq.(28) we see that 
these weak bars end at corotation, as already clear 
from eq.(21), the transition between the oscillatory 
and the unstable regimes coincides with the coro- 
tation condition. The interpretation provided here 
to eq.(22), highlighting the role of orbital diffusion, 
might offer a simple explanation to the commonly 
found truncation of bars at corotation in both or- 
bital structure studies and N-body simulations, e.g. 
the unstable modes of Jalali & Hunter (2005), the 
simulated bars of O'NeiU & Dubinski (2003) or the 
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z=0.05 



z=0.10 




Fig. 2. Disk particles distribution at T = 0. The upper 
panel is 12 x 12 LU, and the other ones are 12 x 4 LU. 

outer of the double bars of Maciejewski & Athanas- 
soula (2008) and Shen & Debattista (2008). 

In summary, weak bars can develop in the os- 
cillatory regime of eq.(18), and in the case of solid 
body rotation curves, present an almost constant 
non-winding flp for a large range of radii, for val- 
ues of Q > 2, perhaps related to the 'pulsating bars' 
seen in some numerical experiments, e.g. Shen & 
Debattista (2007). 

4. NUMERICAL MODELS AND SIMULATIONS 
In order to check the predictions of the model 
presented in the more realistic cases of fully self 



consistent 3D galactic models, we have run a large 
grid of simulations of isolated, stable disk mod- 
els. The models contain three self-consistent compo- 
nents, corresponding to a disk, a bulge and a halo. 
The construction of the models and their following 
evolution were performed using the package Nemo 
(http://carma. astro, umd . edu/ nemo / ) . 

Setting up a stable three-dimensional disk for 
N-body experiments is a difficult task. Sev- 
eral strategies have been suggested including grow- 
ing the disk mass distribution in a self-consistent 
halo/bulge model (Barnes 1988, Athanassoula, Puer- 
ari & Bosma 1997), or treating the spherical compo- 
nents halo-l-bulge as a static background (Sellwood 
& Merritt 1994, Quinn, Hernquist & Fullagar 1993, 
Sellwood 2011). Some authors have also directly 
solved the Jeans equation for the complete disk, 
bulge and halo system to find the velocity disper- 
sions (Hernquist 1993). The more important prob- 
lem in the construction of self-consistent 3 compo- 
nent models is the fragility of a cool disk, generating 
grand design spiral structure and bars, as well as a 
series of other known instabilities (bending modes, 
disk warping, etc). 

Our models were constructed using the Nemo 
routine mkkd95. The routine is based on the method 
described in Kuijken & Dubinski (1995). Their strat- 
egy is the choice of analytic forms for the distribu- 
tion functions (DF) of the three components. For 
the bulge, the DF is taking as a King model (King 
1996). For the halo, the DF is a truncated Evans 
(1993) model for a flattened logarithmic potential. 
For the disk, the more difficult component to build 
up in the model, the authors used a vertical exten- 
sion of the planar DF discussed by Shu (1969) and 
Kuijken & Tremaine (1992). AU the DF's are used 
to calculate the spatial density of each galaxy com- 
ponent, and the Poisson equation 

\ Phalo m 
(28) 

is solved by using spherical harmonic expansion 
by Prendergast & Tomer (1970) with some modifi- 
cations (see details in Kuijken & Dubinski 1995). 

For the evolution of the models, we used the 
Nemo routine gyrf alcON, which is a full-fledged N- 
body code using a force algorithm of complexity 
0{N) (Dehen 2000, 2002), which is about 10 times 
faster than an optimally coded tree code. In the 
program, the gravitational constant is taken equal 
to 1. To match this value, an appropriate normal- 
ization has been chosen with length unit LU = 3 
kpc, and time unit TU = 10^ years. Using these val- 
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12 3 4 



radius 

Fig. 3. Upper panel: tangential velocity curves at 
T = 0. The thick solid line is for the model with 
ch,o = 0.27. The solid, dashed and dotted lines are for 
cr,o = 0.37, 0.47, and 0.57, respectively. Bottom panel: 
Toornre's Q. The lines are as in the upper panel, i.e., 
the thick solid line represents the cooler model, and the 
dotted line represents the warmer one. 

ues, the units of mass, velocity and volume density 
are 6 x 10^° Mq, 293 km s'^, and 2.22 Mq pc'^, re- 
spectively, making the unit of length 3 kpc. In this 
section, we plot all quantities in computational units. 
With gyrf alcON, all the simulations were run using 
a tolerance parameter 9 = 0.5, a softening parameter 
e = 0.05 and a time step of 1/2 These values ensure 
a total energy and angular momentum conservation 
of the order of 10~^ and 10~^, respectively. 

Our models were chosen to be Milky Way 
A (MW-A) like-models from Kuijken & Dubinski 
(1995), having masses of Md ■ Mb : Mh = 0.82 : 
0.42 : 5.2 and number of particles Nd : Nb : Nh = 
160, 000 : 40, 000 : 400, 000. This is twice the parti- 
cle numbers which Kuijken & Dubinski (1995) tested 
in their best model and showed to be in equilib- 
rium. These authors show that disk surface density 
and velocities dispersion profiles hardly change with 



mU.i37.!l)5 P mid_s37_il0 P m=Z.7 
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mM_s37_slO P m=3.7 mkd_s27_zl5 P m=2,7 




13 3 4 5 1 3 3 4 5 

Hadius Radius 

Fig. 4. The P number from eq.(23) as a function of 
radius for 4 of our 16 simulations. Upper panels show the 
models witli aji^a = 0.37 (zq — 0.05, left; zo — 0.1, right), 
and bottom panels show colder models with ajifi = 0.27, 
Zo = 0.1, left; Zo = 0.15, right. The plotted curves are 
the mean P from T = to 200, 2Gyr. The lower curve 
in each panel is for m = 2 and the upper one for m — 7. 
The horizontal lines show limits for P between 1.0 and 
1.15. The upper and lower cuts of these limits and each 
m curve constrain the radial region outwards of which 
we expect the power of that given m to arise. 

time when using the fully self-consistent model with 

300,000 particles; only the disk scale length changes 
of the order of 15%. All our models start with no 
vertex deviation, and result in average values close to 
cee/c/f/f = 0.75 and (Tzz/o'rr = 0.6 for the velocity 
ellipsoid. The fiducial model of Kuijken & Dubinski 
(1995) has a disk central radial velocity dispersion 
Cfl,o = 0.47 and a vertical scale length zq = 0.10. 
We have created a grid of 16 models using aufl = 
0.27, 0.37, 0.47, 0.57 and zq = 0.05, 0.1, 0.15, 0.2, and 
therefore covering the space parameter from cool 
thin disks, to warmer/thicker ones. Here we present 
results for colder /thinner models; in warmer /thicker 
models, the structures present much lower ampli- 
tude, and they are much more chaotic. Our fully self- 
consistent colder /thinner models support low ampli- 
fication of several modes, unlike the unstable models. 
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Fig. 5. Power of the m components as a function of 
radius for the same simulations in figure (Q. Different 
lines are the means for different time intervals: solid line, 
< T < 50, dashed 50 < T < 100, dot-dashed 100 < 
T < 150, and dotted, 150 < T < 200. The hashed areas 
represent the radial range outwards of which the power 
of a given component is expected to increase (1 < P < 
1.15), for each m. For clarity, we added a constant to 
each subsequent m curve above m = 2. 

in which bars and m=2 grand design spiral patterns 
normally form. 

In figure we present a face-on view and edge- 
on snapshots for models with different vertical scale 
length zq. The axisymmetrical aspect of the disk 
remains for all times, only changing about 15% in 
vertical scale. In figure ([3]) we show the tangen- 
tial velocity for models with different disk central 
radial velocity dispersion cr_R,0: a-s well as the initial 
Toomre's Q profiles. 

In figure ^ we show our profiles of P for 4 mod- 
els, 2 for aRfi = 0.37 models (upper panels) and 2 for 
the colder model with cr/j o = 0.27 (bottom panels). 



For each m, the P profile scales with a{R)p{R)~^^^ , 
where cr(R) is the disk radial velocity dispersion pro- 
file and p{R), the disk volume density profile, c.f. 
eq.(23). The plotted curves are the azimuthally av- 
eraged time means of P as a function of radius for 
the total time interval over which the models were 
evolved, from T = to T = 200 = 2Gyr. We 
have draw two horizontal limits for P, P — I (lower 
straight line), and P ~ 1.15 (upper straight line). 
These two limits are used to estimate the radial limit 
outwards of which a given m component is expected 
to show a clear increase in its amplitude, as predicted 
by the developments of the previous sections. 

Now, to calculate the amplitude (or power) of 
each m component as a function of radius, we di- 
vide the disk in 50 rings (from i? = to i? = 5, 
AR = 0.1). For each ring, the ID Fourier transform 
is calculated as: 

1 ^ 

1=1 

where N is the number of disk particles in the 
ring, m is the m-armed component and 9 the polar 
coordinate of each particle in the ring. The real and 
imaginary parts for each m component are 

N N 

Ic(m) — cos{m0i) IsiiTi) = sin{m9i) 

i=l i=l 

(30) 

and the power is then 

Pow^^h{m)^ + Is{m)^ (31) 
The phase of each component is then: 

$m — atan{Is{m) / Ic{m)) /m. (32) 

We have calculated the power and phase of the 
first ten modes, as a function of radius and time. In 
figure (5), we show the power of the to = 2 to 7 
modes, as a function of radius. We present the mean 
value of the power for 4 time intervals. Superim- 
posed on the curves, we have draw shaded areas rep- 
resenting the radial interval resulting from the ana- 
lytic developments of section (2), as inferred through 
eq.(23) applied to figure (4). There is a clear agree- 
ment between the radial positions from figure (U) and 
the radius at which we see an increase of power for 
the different m's, lending empirical support to the 
theoretical model presented in this study. 

Next, we present some plots of phase as a func- 
tion of radius and time in figure (|6]). We present 
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Fig. 6. The phase for the m = 3 (upper panels) and m — 6 (bottom panels) modes for two models appearing in Fig. 
[3 as a function of time and radius (mkd_s27_zl0, left; mkd_s37_z05, right). The phases range from (black) to In/m 
(white). As all features are almost vertical, we see that the perturbations appearing are quite radial, far from tightly 
wound. The horizontal line for each m of each model is the mean of the radial range show in figure 0. Note the clear 
change in the behavior of the phase on crossing the radial position marked, where we expect an increase of power in 
terms of the theory presented. 



the phase of m = 3 (upper panels) and m — 6 (bot- 
tom panels) for the models mkd_s27_zl0 (left) and 
mkd_s37_z05 (right). In these graphs, the phases 
range from to 27r/m. The amplitudes of the per- 
turbations in our isolated models are quite low com- 
pared to those associated with bars and strong m — 2 
spiral arms in unstable models (Puerari, in prepara- 
tion). Even with these low amplitudes, the phase 
shows a clear behaviour, quite ordered, represent- 
ing quasi radial perturbations on the disk, as evi- 
dent from the almost vertical form of the features 
appearing. We also mark in figure (|6]) the mean ra- 
dial position of the corresponding P < 1 threshold 
of eq.(23), as taken from figure At these radial 



positions we notice a change in phase behavior, from 
a more chaotic one (smaller radii) to more ordered 
features (larger radii) . The radial position where this 
transition occurs agrees very well with the theoret- 
ical predictions, with this transition radii occurring 
first for the lower m modes, and only appearing at 
larger radii for higher m modes. 

We end this section with figure (7), which shows 
the regions in flp, R space where power appears for 
the different modes in one of the simulations run. 
We see again a clear agreement with the theoreti- 
cal expectations developed, with the amplitude of 
the modes constrained to the region where ftp < 
rio ± k/to. Comparing with the analytical predic- 
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Fig. 7. The contour plots give the zones where power is 
present for the first 6 modes for one of the models, as a 
function of frequency and radius. The three curves give 
flo and f2o ± K/m 

tions of eq. (21), and given the values of Q of the 
simulated galaxies of ~ 2, we see again good a agree- 
ment. If one wanted to calibrate the details of the 
theory through a comparison of eq.(21) with the re- 
sults of figure(7), we see than a small numerical ad- 
justment, e.g. in the first order scale height estimate 
of eq.(14) of a factor of 2, would furnish an accurate 
accordance of the limit flp predictions of eq.(21) and 
figure (7). 

In summary, we see that even though in the the- 



oretical developments of sections (2) and (3) we ig- 
nored the vertical structure of the galactic disk, the 
radial motion of the stars, and the radical variations 
in surface density, epicycle frequency and velocity 
dispersion, the predictions for the critical radii at 
which quasi-radial perturbations appear, the way in 
which these increase progressively for higher modes, 
and the limit frequencies which these exhibit, are 
all in good accordance with the results of fully self- 
consistent 3D numerical models. 

5. CONCLUSIONS 

We have developed a simple description of a 
quasi-radial density wave in a galactic disk, assum- 
ing only tangential forces due to the perturbation, 
we reach a description of the problem which is com- 
plementary to the classical tight- winding approxima- 
tion. The main features are the substitution of the 
differential rotation for the velocity dispersion of the 
disk in the resulting dispersion relation. The result- 
ing dispersion relation allows a clear understanding 
of why lower m modes are more unstable than higher 
m modes, and of why the higher m values appear 
progressively at larger radii. 

Weak bars in harmonic potential central regions 
appear in the oscillatory regime, an expression for flp 
is found which shows a divergence towards R = 0, 
perhaps helping to understand the phenomenon of 
increasingly rapidly rotating bars within bars seen 
in numerical simulations. The extent of this weak 
bars is seen to be naturally limited by the corotation 
condition. 

The expectations for quasi-radial spiral arms are 
tested in detail through extensive numerical simu- 
lations, which confirm to good accuracy the main 
results of the theory developed here. 
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